par(mfrow=c(2, 2))
data1_15 <- read.csv("data/时间序列分析——基于R（第2版）案例数据/csv/A1_15.csv")
ts_fertility <- ts(data1_15$fertility)

# 时序图
plot(ts_fertility, type="o", pch=5, col='#39CBB4')


# ADF检验
# install.packages("aTSA")
library(aTSA)
adf.test(ts_fertility, nlag=3)

# 一阶差分
dif_fertility<-diff(ts_fertility)
plot(dif_fertility, type="o", pch=5, col='#39CBB4')

# ADF检验
adf.test(dif_fertility, nlag=3)


# 白噪声，纯随机检验
for( k in 1:3) print(Box.test(dif_fertility, lag=6*k, type="Ljung-Box"))


# 自相关图
acf(dif_fertility, lag.max=20)
pacf(dif_fertility, lag.max=20)

# ARIMA((1, 4), 1, (1, 4, 5))
model11 <- arima(ts_fertility, order=c(4, 1, 0), method="ML", transform.pars=F, fixed = c(NA, 0, 0, NA))
model12 <- arima(ts_fertility, order=c(0, 1, 5), method="ML", transform.pars=F, fixed = c(NA, 0, 0, NA, NA))
model11
model12

# 模型显著性检验
ts.diag(model11)
ts.diag(model12)
# 参数显著性检验
# 粗略的
model11
model12
# 精确的 构造t统计量，求P值
t = abs(model11$coef)/sqrt(diag(model11$var.coef))
pt(t, length(ts_fertility)-length(model11$coef), lower.tail = F)

t = abs(model12$coef)/sqrt(diag(model12$var.coef))
pt(t, length(ts_fertility)-length(model12$coef), lower.tail = F)


# 或者
# install.packages("forecast")
library(forecast)
model21 <- Arima(ts_fertility, order=c(4, 1, 0), method="ML", include.drift = T, fixed = c(NA, 0, 0, NA))
model22 <- Arima(ts_fertility, order=c(0, 1, 5), method="ML", include.drift = T, fixed = c(NA, 0, 0, NA, NA))
model21
model22
